function [ VAR ] = rfvar(VAR)
%rfvar Estimate RF var parameters.
%  Input:
%   Y = T x Ny matrix of LHS variables
%   X = T x (Nx-1), matrix of RHS variables excluding constant  
%       X(t)=[y(t-1)',y(t-2),...,y(t-p)]; 
%  Output:
% bet = Nx x Ny matrix of VAR OLS parameter estimates ("beta")
%   S = Ny x Ny matrix of sum of squared residuals
% Sigma= S/T covariance matrix of RF residuals
% betvec = Nx*Ny x 1 vectorized matrix of bet

% Estimate PHI, S
Y   = VAR.Y; 
X   = VAR.X;
T   = size(Y,1);
Ny  = size(Y,2);
X1  = [X ones(T,1)]; % add constant column to X
Nx  = size(X1,2);
bet = (X1'*X1)\(X1'* Y);
res = Y - X1*bet;
S   = res'*res;
Sigma=S./(T - Nx);
betvec= reshape( bet, [Nx*Ny , 1]);  % vectorize estimates

% store in "VAR" structure
VAR.T    = T;
VAR.Ny   = Ny;
VAR.Nx   = Nx;
VAR.bet  = bet;
VAR.res  = res;
VAR.S    = S;
VAR.Sigma= Sigma;
VAR.betvec= betvec;
end

